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ABSTRACT 

With increasingly large data sets, weak lensing measurements are able to measure 
cosmological parameters with ever greater precision. However this increased accuracy 
also places greater demands on the statistical tools used to extract the available in- 
formation. To date, the majority of lensing analyses use the two point-statistics of the 
cosmic shear field. These can either be studied directly using the two-point correlation 
function, or in Fourier space, using the power spectrum. But analyzing weak lensing 
data inevitably involves the masking out of regions for example to remove bright stars 
from the field. Masking out the stars is common practice but the gaps in the data 
need proper handling. In this paper, we show how an inpainting technique allows us 
to properly fill in these gaps with only A/" log A/" operations, leading to a new image 
from which we can compute straight forwardly and with a very good accuracy both 
the pow er spectrum and the bispectrum. We propose then a new method to compute 
the bispectrum with a polar FFT algorithm, which has the main advantage of avoiding 
any interpolation in the Fourier domain. Finally we propose a new method for dark 
matter mass map reconstruction from shear observations which integrates this new 
inpainting concept. A range of examples based on 3D N-body simulations illustrates 
the results. 
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1 INTRODUCTION 

The distortion of the images of distant galaxies by gravita- 
tional lensing offers a direct way of probing the statistical 
properties of the dark matter distribution in the Universe; 
without making any assumption about the relation between 
dark a nd visible m at ter, see iBartelmann and Schneider 
(200l');'Mellier' ("iQQgl lJVan Waerbeke etHI (|200lh : iMelliCT 
(2002); Refregier (2003|). This weak lensing effect has been 
detected by several groups to derive constraints on cosmo- 
logical parameters. Analyzing an image for weak lensing in- 
volves inevitably the masking out of regions to remove bright 
stars from the field. Masking out the stars is common prac- 
tice but the gaps in the data need proper handling. 

At present, the majority of lensing analyses use the two 
point-statistics of the cosmic shear field. These can either 
be st udied directly using the two-point correlation func- 



tion (iMaoli et al.1 boOl]: [ Refregier et al.1 l2002l : iBacon et all 
I2OO3I : iMassev et al.1 l2005ir or in Fourier space, using the 



power spectrum ( Brown et al.|[2003l V Higher order statisti- 
cal measures, such as three or four-point correlation func- 



tions have been st udied ("B ernardeau et al" I l2003l :l Pen et al. 
I2OO3I : IJarvis et al.1 12004) and have shown to provide addi- 
tional constraints on cosmological parameters. 

Direct measurement of the correlation function, through 
pair counting, is widely used since this method is not bi- 
ased by missing data, for instance the ones arising from the 
masking of bright stars. However, this method is computa- 
tionally intensive, requiring 0{N^) operations. It is there- 
fore not feasible to use it for future ultra-wide lensing sur- 
veys. Measuring the power spectrum is significantly less de- 
manding computationally, requiring 0(A^log A^) operations, 
but is strongly affected by missing data. The estimation of 
power spectra from various types of data is becoming in- 
creasingly important also in other cosmological applications 
such as in the analysis of Cosmic Microwave Background 
(CMB) temperature maps or galaxy clustering. In the liter- 
ature, a large number of papers have appeared in the last few 
years that discuss various solutions to the problem of power 
spectrum esti mation from large data sets with complete or 
missing data (' Bond et al.lll998l : iRuhl et al.ll2003l : iTegmarkI 
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19971: iHivon et al J '2002"; 'Ha nsen et aD l2002l; ISzapudi et aD 
2001bl lal: [Efetathi ou 2004; Sza pudi et aLll2005l ). As explained 
in details in section I3.4[ they present however some limi- 
tations such as numerical instabilities which require to to 
regularize the solution. In this paper, we investigate an al- 
ternative approach based recent work in harmonic analysis 



A new approach: inpainting 

Inpainting techniques are well known in the image process- 
ing litterature and are used to fill the gaps (i.e. to fill the 
missing data) by inferring a maximum information from the 
remaining data. In other words, it is an extrapolation of the 
missing information using some priors on the solution. We 
investigate in this paper how to fill-in judiciously masked 
regions so as to reduce the impact of missing data on the 
estimation of the power spectrum and of higher order sta- 
tistical measures. The inpainting approach we propose re- 
lies on a long-standing discipline in s tatistical estimation 
the ory; estimatio n with missing data (Dempster et al. I ll977l : 
[Little and Rubin! 1 19871 ). We propose to use an inpainting 
method that r elies on the spars e representation of the data 
introduced bv'E lad et al.l (|2005l l. In this work, inpainting is 
stated as a linear inverse ill-posed problem, that is solved 
in a principled bayesian framework, and for which the pop- 
ular Expectation-Maximization mechanism comes as a nat- 
ural iterative algo rithm because of physically missing data 
(jFadiU et al. l l2007l ). Doing so, our algorithm exhibits the fol- 
lowing advantages: it is fast, it allows to estimate any statis- 
tics of any order, the geometry of the mask does not imply 
any instability, the complexity of the algorithm does not de- 
pend on the mask nor on data weighting. We show that, 
for two different kinds of realistic mask (similar to that for 
CFHT and Subaru weak lensing analyses), we can reach an 
accuracy of about 1% and 0.3% on the power spectrum, and 
an accuracy of about 3% and 1% on the equilateral bis- 
pectrum. In addition, our method naturally handles more 
complicated inverse problems such as the estimation of the 
convergence map from masked shear maps. 

This paper is organized as follows. Section [2] describes 
the simulated data that will be used to validate the pro- 
posed methods, especially how large statistical samples of 
3D N-body simulations of density distribution have been 
produced using a grid architecture and how 2D weak lensing 
mass maps have been derived. It also shows typical kinds of 
masks that need to be considered when analysing real data. 
Section [3] is an introduction to different statistics which are 
of interest in weak lensing data analysis and a fast and accu- 
rate algorithm to compute the equilateral bispectrum using 
a polar Fast Fourier Transform (FFT) is introduced. The 
speed of the bispectrum algorithm arises from the quickness 
of the polar Fourier transform and the accuracy comes from 
the output polar grid of the polar Fourier transform, thus 
avoiding the interpolation of coefficients in Fourier space. In 
section m we present our inpainting method for gap filling in 
weak lensing data, and we show that it is a fast and accurate 
solution to the missing data problem for second and third 
order statistics calculation. In section O we propose a new 
approach to derive dark matter mass maps from incomplete 
weak lensing shear maps which uses the inpainting concept. 
Our conclusions are summarized in section [6] 



2 SIMULATIONS OF WEAK LENSING MASS 
MAPS 

2.1 3D N-body cosmological simulations 

We have run realistic simulated convergence mass maps de- 
rived from N-body cosmological simulations using the RAM- 
SES code (Teyssier 2002). The cosmological model is taken 
to be in concordance with the ACDM model. We have cho- 
sen a model with the following parameters close to WMAPl 
: Qrn — 0.3, (78 — 0.9, = 0.7, h — 0.7 and we have run 
33 realizations of the same model. Each simulation has 256^ 
particles with a box size of 162 Mpc/h. We refined the base 
grid of 256^ cells when the local particle number exceeds 10. 
We further refined similarly each additional levels up to a 
maximum level of refinement of 6, corresponding to a spa- 
tial resolution of 10 kpc/h, making certain the particle shot 
noise remains at an acceptable level (see Fig. [14] and Fig. 

[El). 

2.2 Grid computing 

The simulation suite was deployed on a gr id architecture de - 
signed under the project name Grid' 5000 ( Bolze et al.ll2006l l. 
For that purp ose, we use a newly de veloped middle-ware 
called DIET ( Caron and Despred 12006. ) in order to compile 
and execute the RAMSES code on a widely inhoraogen eous 
grid system. For this experiment (jCaniou et al.1 l2007l ). 12 
clusters have been used on 7 sites for a duration time of 48 
hours. In total, 816 grid nodes have been used for the present 
experiment, leading to the execution of 33 complete simu- 
lations. Note that this overall computation was completed 
within 2 days. This would have taken more than one month 
on a regular 32 processor cluster. 

2.3 2D Weak lensing mass map 

In N-body simulations, which are commonly used in cos- 
mology, the dark matter distribution is represented using 
discrete massive particles. The simplest way to deal with 
these particles is to map their positions onto a pixelised 
grid. In the case of multiple sheet weak lensing, we do this 
by taking slices through the 3D simulations. These slices 
are then projected into 2D mass sheets. The effective con- 
vergence can subsequently be calculated by stacking a set of 
these 2D mass sheets along the line of sight, using the lens- 
ing efficiency f unction. This is a proc edure that has been 
used before bv IVale and White! (|2003h . where the effective 
2D mass distribution K,e is calculated by integrating the den- 
sity fluctuation along the line of sight. Using the Born ap- 
proximation which neglects the facts that the light rays do 
not follow straight lines, the convergence (see eq. [16]) can be 
numerically expressed by : 

^ ^ Xoa(xO I " ^''^) 

where Hq is the Hubble constant, Qrn is the density of mat- 
ter, c is the speed of light, L is the length of the box x ^tre 
co-moving distances, with xo being the co-moving distance 
to the source galaxies. The summation is performed over the 
i^^ box. The number of particles associated with a pixel of 
the simulation is rip, the total number of particles within a 
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Figure 1. Simulated weak lensing convergence map for the 
ACDM cosmological model with ag = 0.9 and = 0.3. The 
region shown is 1° x 1°. 




Figure 2. Left, mask pattern of Subaru Survey 0.575° x 0.426° 
(with SuprimeCam camera) right, mask pattern of CFHTLS data 
on Dl field 1° x 1° (with the MegaCam camera). Courtesy Joel 
Berge. 



simulation is Nt and s — Lp/L where Lp is the length of the 
plane doing the lensing. R is the size of the 2D maps and 
Ar/i = Q^YLd r2 are co-moving distances. 

We have derived 100 simulated weak lensing mass maps 
from the previous 3D simulations. Fig. [1] shows a zoom of 
one of the 2D maps obtained by integration of the 3D den- 
sity fluctuation on one of these realizations. The total field 
is 1.975 X 1.975 square degrees, with 512 x 512 pixels and 
we assume that the sources lie at exactly z = 1. The over- 
densities (peaks) correspond to halos of groups and clusters 
of galaxies. The typical standard deviation values of k are 
thus of the order of a few percent. 



2.4 2D Weak lensing mass map with missing data 

Loss of data can be caused by many factors. Missing data can 
be due to camera CCD defect or from bright stars in the field 
of view that saturate the image around them as seen in weak 
lensing surveys with diff e rent telescope s (Hoekstr a et al 



iensmg surveys witn din e rent telescope s ( MoeKstr a et al 
20061 : iHamana etHI l2003l : iMassev et aP [2005 : Be rge et al 



2008 



Fig [2] shows the mask pattern of CF HTLS image in the 
Dl-field with about 20 % of missing data (iBerge et al.ll2008h 
and that of SUBARU image covering a part of the same 
field with about 10 % of missing data. The mask pattern 
depends essentially on the field of view and on the quality 
of the optics. 

In our simulations, we have chosen to consider these two 
typical mask patterns to study the impact of gaps in weak 
lensing analysis. The problem is now to extract statistical 
informations from weak lensing data with such masks. 



3 STATISTICS 

3.1 Two-point statistics 

Statistical weak gravitational lensing on large scales probes 
the projected density field of the matter in the Universe : 
the convergence k,. Two-point statistics have become a stan- 
dard way of quantifying the clustering of this weak lensing 
convergence field. Much of the interest in this type of anal- 
ysis comes from its potential to constrain the spectrum of 
density fluctuations present in the late Universe. All second- 
order statistics of the convergence can be expressed as func- 
tions of the two-point correlation function of k or its Fourier 
transform, the Power Spectrum P^. 

• Two-point correlation function : 
Direct two-point correlation function estimators are 
based on the notion of pair counting. As a result of the 
Universe's statistical anisotropy, it only depends on |^| the 



distance between the position 
given by : 



and the position Oj , and is 



N N 

i=l i=i 



(2) 



where No is the number of pairs separated by a distance 
of 1^1 . Its naive implementation requires 0(A/'^) operations. 
Pair counting can be sped up, if we are interested in mea- 
suring the correlation function only on small sca l es. In that 
case the double-tree algorithm by iMoore et al.l (|200lh re- 
quires approximatively 0{N log N) operations. However if 
all scales are considered, the tree-based algorithm slows 
down to 0(N'^) operations like naive counting. 

• Power spectrum : 
The power spectrum P^^ is the Fourier transform of the two- 
point correlation function (by the Wiener- Khinchine theo- 
rem). Because of the rotational invariance derived from the 
Universe isotropy, the Fourier transform becomes a Hankel 
transform : 

p4q) = ^ j^^ c{e)M2nqe)ede, (s) 

where Jo is the zero order Bessel function. Computationally, 
we can estimate the power spectrum directly from the signal: 



(4) 



where k denotes Fourier transform of the convergence. Thus 
we can take advantage of the FFT algorithm to quickly es- 
timate the power spectrum. 

• Sensitivity to missing data: 
In weak lensing data analysis, it is common practice to mask 
out bright stars, which saturate the detector. This requires 
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an appropriate post-treatment of the gaps. Contrarily to the 
two-point correlation function, the power spectrum estima- 
tion is strongly sensitive to missing data. Gaps generate a 
loss of power and gap edges produce distortions in the spec- 
trum that depend on the size and the shape of the gaps. 



3.2 Three-point statistics 

Analogously with two-point statistics, third-order statistics 
are related to the three-point correlation function of k or its 
Fourier transform the Bispectrum B^. It is well established 
that the primordial density fluctuations are near Gaussian. 
Thus, the power spectrum alone contains all information 
about the large-scale structures in the linear regime. How- 
ever, gravitational clustering is a non-linear process and in 
particular, on small scales, the mass distribution is highly 
non-gaussian. Three-point statistics are the lowest-order 
statistics to quantify non-gaussianity in the weak lensing 
field and thus provides additional information on structure 
formation models. 

• Three-point correlation function : 
Direct three-point correlation function estimators C^^k are 
based on the notion of triangle counting. It depends on dis- 
tances di, d2 and ds between the three spatial positions Oi, 
Oj and Ok of the triangle vertices formed by three galaxies, 
and is given by : 



(5) 



where < . > stands for the expected value. The naive imple- 
mentation requires 0{N^) operations and can consequently 
not be considered on future large data sets. One configu- 
ration, that is often used, is the equilateral configuration, 
wherein di = d2 = ds = d. The three-point correlation func- 
tion can then be plotted as a function of d. The configuration 
dependence being weak (|Coorav and Huir200lh , the equilat- 
eral configuration has become standard in weak lensing : 
first, because of its direct interpretation and second because 
its implementation is faster. The equilateral three-point cor- 
relation function estimation can be written as follows : 



1 



N N N 

i=l j=l k=l 



(6) 



where Nd is the number of equilateral triangles whose side is 
d. Whereas primary three-point correlation implementation 
requires 0{N^) operations, equilateral triangle counting can 
be sped up to 0{N^) operations. But this remains too slow 
to be used on future large data sets. 

• Bispectrum : 
The complex Bispectrum is formally defined as the Fourier 
transform of the third-order correlation function. We assume 
the field k to be statistically isotropic, thus its bispectrum 
only depends on distances |/c2| and Ifel: 



B{\kil \k2l |fe|) oc < ^(|/ci|)^(|/c2|)r (Ifel) > . 



(7) 



If we consider the standard equilateral configuration, the tri- 
angles have to verify : ki = k2 = ks = k and the bispectrum 
only depends on k : 



Biky"" oc < k{k)k{k)k*(k)> 

• Sensitivity to missing data : 



(8) 




Figure 3. Equilateral bispectrum configuration in Fourier Space. 
Equilateral triangles must be inscribed in a circle of origin (0, 0) 
and of radius k. 



Like two-point statistics, the three-point correlation func- 
tion is not biased by missing data. On the contrary, the 
estimation of bispectrum is strongly sensitive to the missing 
data that produce important distortions in the bispectrum. 
For the time being, no correction has been proposed to deal 
with this missing data on bispectrum estimation and the 
three-point correlation function is computationally too slow 
to be used on future large data sets. 



3.3 The weak lensing statistics calculation from 
the polar FFT 

Assuming the gaps are correctly filled, the field becomes 
stationary and the Fourier modes uncorrelated. We present 
here a new method to calculate the power spectrum and the 
equilateral bispectrum accurately and efficiently. 

To compute the bispectrum, we have to average over 
the equilateral triangles of length k . Fig. [3] shows the form 
that such triangles in Fourier space. For each /c, we inte- 
grate over all the equilateral triangles inscribed in the circle 
of origin (0,0) and of radius \k\. Because of the rotational 
symmetry we have only to scan orientation angles from to 
27r/3. The bispectrum is obtained by multiplying the Fourier 
coefficients located at the three vertices. 

Similarly, to compute the power spectrum, we have to 
average the modulus squared of the Fourier coefficients lo- 
cated in a circle of origin (0, 0) and of radius k. 



The polar Fast Fourier transform (polar FFT) 

It requires some approximations to interpolate the Fourier 
coefficients in an equi-spaced Cartesian grid, as shown in 
Fig. |3]on the left. In order to avoid these approximations, a 
solution consists in using a recent method, called polar Fast 
Fourier Transform that is a powerful tool to manipulate the 
Fourier transform in polar coordinates. 

A fast and ac c urate Polar FFT has been proposed 
by lAverbuch et J] diooEh . For a given two-dimensional 
signal of size N, the proposed algorithm's complexity is 
O(iVlogiV), just like in a Cartesian 2D-FFT. The polar 
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Figure 4. Calculation of the mean power per frequency using a 
regular grid (left) and using a polar grid (right) 



FFT is just a particular case of the more general prob- 
lem of finding the Fourier transform in a non-equispaced 
grid (|Keiner et all l2006l ). We have used the NFFT (Non 
equi-spaced Fast Fourier Transform, software available at 
http:/ /www-US er.tu-chemnitz.de/^ potts /nfft) to compute a 
very accurate power spectrum and equilateral bispectrum. 
Fig. m (right) shows the grid that we have chosen. For each 
radius, we have the same number of points. The calculation 
of the average power associated to each equilateral triangle 
along the radius becomes easy and approximations are no 
longer needed. 



Algorithms 

The Polar-FFT bispectrum algorithm is: 



1. Forward polar Fourier transform of the convergence hi. 

2. Set the radius (in polar coordinates) r to 0. Iterate: 
3. Set the angle (in polar coordinates) 6 to 0. Iterate: 

4. Locate the cyclic equilateral triangle whose one vertex 
(or corner) have (r, 0) as coordinate. A cyclic triangle is 

a triangle inscribed in a circle it means the sides are chords 
of the circle. 

5. Perform the product of the Fourier coefficients located 
at the three corners of the cyclic equilateral triangle. 

6. 6* = 6> + (56* and if ^ < 27r/3 return to step 4. 

7. Average the product over all the cyclic equilateral triangle 
inscribed in the circle of radius r. 

8. r = r -\- 1 and if r < rmax, return to Step 3. 



The Fourier coefficients at the three corners of the 
cyclic equilateral triangle are easy to locate using a polar 
grid. We don't need to interpolate to obtain the Fourier 
coefficient values as we have to do with a Cartesian grid. 
In addition to be accurate, this computation is also very 
fast. Indeed, in the simulated field that covers a region 
of 1.975° X 1.975° with 512 x 512 pixels, using a 2.5 
GHz processor PC-linux, about 60 seconds are needed to 
complete the calculation of the equilateral bispectrum and 
the process only requires 0{N\ogN) operations. The bulk 
of computation is invested in the polar FFT calculation. 
This algorithm will be used in the experiments in section ^ 

Similarly the Polar-FFT power spectrum algorithm is: 



1. Forward polar Fourier transform of the convergence k,. 

2. Take the modulus squared of the polar Fourier transform of 
the convergence. 

3. Set the radius (in polar coordinates) r to 0. Iterate: 

4. Average the power over all the possible angles in the circle of 
radius r. 

5. r = r -\- 1 and if r < rmax, return to Step 4. 



3.4 The missing data problem: state of the art 

Second Order Statisitics 

In the literature, a large number of studies have been pre- 
sented that discuss various solutions to the problem of power 
spectrum estimation from large data sets with complete or 
missing data. These papers can be roughly grouped as fol- 
lows: 

• Maximum likelihood (ML) estimator: two types of ML 
estimators have been discussed. The first uses an iterative 
Newton-type algorithm to maximize the likelihood score 
function without any assumed explicit form on the co- 
variance matrix of the observed data ([Bond et al.l Il998l : 
iRuhl et aDbOOSi ) . The sec ond one is based on a model of the 
power- spectrum; see e.g. lTegmarl3 (|l997f ). The ML frame- 
work allows to claim optimality in ML sense and to derive 
Cramer-Rao lower-bounds on the power spectrum estimate. 
However, ML estimators can become quickly computation- 
ally prohibitive for current large-scale datasets. Moreover, to 
correct for masked data, ML estimators involve a "deconvo- 
lution" (analogue to the PPS method below) that requires 
the inversion of an estimate of the Fisher information ma- 
trix. The latter depends on the mask and may be singular 
(semidefinite positive) as is the case for large galactic cuts, 
and a regularization may need to be applied. 

• P seudo power-spectrum (PPS) estimators: iHivon et al.l 
(|2002l ) proposed the MASTER method for estimating 
power-spectra for both Cartesian and spherical grids. Their 
algorithm was designed to handle missing data such as the 
gala ctic cut in CMB dat a using apodization windows; see 
also iHansen et aD (|2002 h. These estimators can be evalu- 
ated efficiently using fast transforms such as the spherical 
harmonic transform for spherical data. Besides their fast im- 
plementation, PPS-based methods also allow us to derive an 
analytic covariance matrix of the power spectrum estimate 
under certain simplifying assumptions (e.g. diagonal noise 
matrix, symmetric beams, etc.). However, the deconvolu- 
tion step in MASTER requires the invertion of a coupling 
matrix which depends on the power spectrum of the apodiz- 
ing window. The singularity of this matrix strongly relies on 
the size and the shape of missing areas. Thus, for many 
mask geometries, the coupling matrix is likely to become 
singular, hence making the deconvolution step instable. To 
cope with this limitation, one may reso rt to regulari z ed in - 
verses. This is for instance the case in iHivon et al.1 (|2002l ) 
where it is proposed to bin the CMB power spectrum. Do- 
ing so, the authors implicitly assume that the underlying 
power spectrum is piece- wise constant, which yields a loss 
of smoothness and resolution. 

A related class of sub-optimal estimators use fast evalu- 
ation of the two-point correlation function, which can then 
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be transformed to give an estimate of the power spectrum. 
M ethods of this t y pe (use d e.^. fo r the CMB) are described 
by Szapudi et all ((2001bl fl l200 5h (the SPICE method and 
its Euchdean version eSPICE). This class of estimators is 
closely related, though not exactly equivalent to the PPS 
estimator. However, t here are two is s ues wi t h the estima- 
tion formulae given bv lSzapudi et al.1 (|2001bl . I2OO5I I: 

The first one concerns statistics. SPICE uses the Wiener- 
Khinchine theorem in order to compute the 2PCF in the 
direct space by a simple division between the inverse Fourier 
transform of the (masked) data power spectrum and the 
inverse Fourier transform of the mask power spectrum. But 
when data are masked or apodized, the resulting process is 
no longer wide-sense stationary and the Wiener-Khinchine 
theorem is not strictly valid anymore. 

The second one is methodological. Indeed, to correct for 
missing data, instead of inverting the coupling matrix in 
spherical harmonic or Four ier spaces as done in MASTER 
(the " deconvolution step" ) , ISzapudi et al.l (|2005h suggest to 
invert the coupling matrix in pixel space. They accomplish 
this by dividing the estimated auto-correlation function of 
the raw data by that of the mask. But, this inversion (decon- 
volution) is a typical ill-posed inverse problem, and a direct 
division is unstable in general. This could be alleviated with 
a regularization scheme, which needs then to be specified for 
a given application. 

MASTER has been designed for data on the sphere, and 
no code is available for Cartesian maps. eSPICE has been 
proposed for computing the power spectrum of a set points 
(e.g. galaxies), and has not tested for maps where each pixel 
position has an associated value (i.e. weight). Therefore a 
pubUc code for cartesian pixel maps remains to be devel- 
oped. 

• Other estimators: some ML estimators that make use 
of the scanning geometry in a specific experiment were pro- 
posed in the literature. A hybrid estimator was also proposed 
that combines an ML estim ator at low multip oles and PPS 
estimate at high multipoles (|Efstathioull2Q04l ). 

The interested reader may refer to lEfstathioul (|2004l l 
for an extended review, further details and a comprehensive 
comparative study of these estimators. 



Third Order Statisitics 

The ML estimators discussed above heavily rely on the 
Gaussianity of the field, for which the second-order statis- 
tics are the natural and sufficient statistics. Therefore, to 
estimate higher order statistics (e.g. test whether the pro- 
cess contains a non-gaussian contribution or not) , the strat- 
egy must be radically changed. Many authors have al- 
ready addressed the problem of three-point statistics. In 
iKilbinger and Schneider (2005), the authors calculate, from 
ACDM ray tracing simulations, third-order aperture mass 
statistics that contain information about the bispectrum. 
Many authors have already derived analytical predictions for 
the three-po i nt corre lati on function and the bispectrum (e.g . 
Ma and Frvl l2000al|3 ; IScoccimarro and CouchmanI I2OOII : 
Coorav and Hul l200lh . Estimating three -point correlation 



function from data has already be done (|Bernardeau et al] 
I2OO2I ) but can not be considered in future large data sets 
because it is computationally too intensive. In the conclu- 



sion of ISzapudi et al.l (|2001bl l. the authors briefly suggested 
to use the p-point correlation functions with implementa- 
tions that are at best 0(iV(logiV)^"^). However, it was not 
clear if this sug gestion is valid for the missing data case. 
IScoccimarro et al.l ((l998l ) proposed an algorithm to com- 
pute the bispectrum from numerical simulations using a Fast 
Fourier transform but without consi dering the case of inco m- 
plete data. This method is used by Fosalba et al. (2005) to 
estimate the bispectrum from numerical simulations in or- 
der to compare it with the analytic halo mo del predictions. 
Some recent stud i es on CMB real data (|Komatsu et al.l 
I2OO5I : lYadav et al.l l2007l l concentrate on the non-gaussian 
quadratic term of primordial fluctuations using a bispectrum 
analysis. But correcting for missing data the full higher- 
order Fourier statistics still remain an outstanding issue. 
In the next section, we propose an alternative approach, 
the inpainting technique, to derive 2nd order and 3rd order 
statistics and possibly higher-order statistics. 



4 WEAK LENSING MASS MAP INPAINTING 

Here we describe a new approach, based on the inpainting 
concept. 



4.1 Introduction 

In order to compute the different statistics described pre- 
viously, we need to correctly take into account the miss- 
ing data. We investigate here a new approach to deal with 
the missing data problem in weak lensing data set, which is 
called inpainting in analogy with the recovery process muse- 
ums experts use for old and deteriorated artwork. Inpainting 
techniques are well known in the image processing litterature 
and consist in filling the gaps (i.e. missing data). In other 
words, it is an extrapolation of the missing information using 
some prior on the sol ution. For instance, Guillermo Sapiro 
and his coll aborators (|Ballester et al.ll200ll : iBertalmio et al.l 
I2OOI", 2000) use a prior relative to a smooth continuation 
of isophotes. This principle leads to nonlinear partial dif- 
ferential equation (PDE) model, propagating information 
from the boundari es of the holes w hile guarant eeing smooth- 
ness i n some way (IChan and Shen 2001: Masn ou~and Morel 
|2002; Bor nemann a nd Mar j I2OO6I : [Chan et al. 2006). Re- 
cently, Elad et al.l f|2005) introduced a novel inpainting al- 
gorithm that is capable of reconstructing both texture and 
smooth image contents. This algorithm is a direct exten- 
sion of the MCA (Morphological Component Analysis), de- 
signed for the separ ation of an image into different semantic 
components (Starck et al.ll2005l : Istarck et al.ll2004h . The ar- 
guments supporting this method were borrowed from the 
theory of Compressed Sensi ng (C S ) rece n tly d evelope d by 
Donoho (2004) and Candes et al.1 ((20041 ): [Candes and Taol 
(|2005l . I2OO4I ). This new method uses a prior of sparsity in 
the solution. It assumes that there exists a dictionary (i.e. 
wavelet. Discrete Cosine Transform, etc) where the complete 
data is sparse and where the incomplete data is less sparse. 
For example, mask borders are not well represented in the 
Fourier domain and create many spurious frequencies, thus 
minimizing the number of frequencies is a way to enforce the 
sparsity in the Fourier dictionary. The solution that is pro- 
posed in this section is to judiciously fill-in masked regions 
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so as to reduce the impact of missing data on the estima- 
tion of the power spectrum and of higher order statistical 
measures. 



4.2 Inpainting based on sparse decomposition 

The classical image inpainting problem can be defined as 
follows. Let X be the ideal complete image, Y the observed 
incomplete image and M the binary mask (i.e. Mi = 1 if we 
have information at pixel i, Mi — otherwise). In short, we 
have: Y — MX. Inpainting consists in recovering X knowing 
Y and M. 

In many applications - such as compression, de-noising, 
source separation and, of course, inpainting - a good and 
efficient signal representation is necessary to improve the 
quality of the processing. All representations are not equally 
interesting and there is a strong a priori for sparse repre- 
sentation because it makes information more concise and 
possibly more interpret able. This means that we seek a rep- 
resentation a = X of the signal X in the dictionary $ 
where most coefficients ai are close to zero, while only a few 
have a significant absolute value. 

Over the past decade, traditional signal representations 
have been replaced by a large number of new multiresolution 
representations. Instead of representing signals as a super- 
position of sinusoids using classical Fourier representation, 
we now have many available alte rnative dictionaries such as 
wavelets (Mallat 1989 ), ridgelets (ICandes and Donohdll999l l 



wavelets (Mallat 19^9 ), ridgelets ( (Jandes and JJonolioll 19991 ) 
or curvelets (Starck et al. I I2OO3I : ICandes et al.ll2006h . most 
of which are overcomplete. This means that some elements 
of the dictionary can be described in terms of other ones, 
therefore a signal decomposition in such a dictionary is not 
unique. Although this can increase the complexity of the sig- 
nal analysis, it gives us the possibility to select among many 
possible representations the one which gives the sparsest rep- 
resentation of our data. 

To find a sparse representation and noting ||z||o the /o 
pseudo-norm, i.e. the number of non-zero entries in z and 
ll^ll the classical I2 norm (i.e. \\z\\^ — ^^^(^fc)^), we want to 
minimize: 



mm 

X 



subject to II y - MX a, (9) 



where a stands for the noise standard deviation in the noisy 
case. Here, we will assume that no noise perturbs the data 
y, cr = (i.e. the constraint becomes an equality). As dis- 
cussed later, extension of the method to deal with noise is 
straighforward. 

It has also been shown that if X is sparse enough, 
the /o pseudo-norm can al so be replaced by the c onvex h 
norm (i.e. \\z\\i = Efc Nfcl) (Donoh o and Hudl200lh . The so- 
lution of such an optimisation task can be obtain ed thr ough 
an it erative thresholding algorithm called MCA (lElad et al.l 
I2OO5I I : 



X 



n+l 



^ A$,A^(x" + M(y-x")), 



(10) 



where the nonlinear operator /\<^^\[Z) consists in: 



• decomposing the signal Z on the dictionary $ to derive 
the coefficients a — 

• threshold the coefficients: a = p(a,A), where the 
thresholding operator p can either be a hard thresholding 
(i.e. p{ai,X) = ai if \ai\ > A and otherwise) or a soft 



thresholding (i.e. p{ai,X) = sign(ai)max(0, \ai\ — A)). The 
hard thresholding corresponds to the lo optimization prob- 
lem while the soft-threshold solves that for h. 

• reconstruct Z from the thresholded coefficients a. 

The threshold parameter An decreases with the iteration 
number and it plays a part similar to the cooling pa- 
rameter of the simulated annealing techniques, i.e. it al- 
lows the solution to escape from local minima. More de- 
tails relati ve to this o ptimizatio n problem can be found in 
Combett es and WaisI jioOS); F adih et al.1 (|2QQ7l ). For many 
dictionaries such as wavelets or Fourier, fast operators ex- 
ist to decompose the signal so that the iteration of eq. [10] 
is fast. It requires us only to perform, at each iteration, a 
forward transform, a thresholding of the coefficients and an 
inverse transform. The case where the dictionary is a union 
of subdictionaries $ = {$1, $t} wh ere each ^i has a fas t 
operator has also been investiga ted in [Starck et~al ] (|2004h : 
Elad et al. (2005); Fadih et al. (2^3)- We will discuss in the 
following the choice of the dictionary for the weak lensing 
inpainting problem. 

In general, there are some restrictions in the use of in- 
painting based on sparsity that arise from the link between 
the sparse representation dictionary and the masking oper- 
ator M. The first restriction is that the proposed inpainting 
method based on sparse representation assumes that a good 
representation for the data is not a good representation of 
the gaps. This means that features of the data need few coef- 
ficients to be represented, but if a gap is inserted in the data 
a large number of coefficients will be necessary to account 
for this gap. Then by minimizing the number of coefficients 
among all the possible coefficients, the initial data can be ap- 
proximated. Secondly, the inpainting is possible if the gaps 
are smaller than the largest dictionary elements. Indeed, if 
a gap removes a part of an object that is well represented 
by one element of the dictionary, this object can be recov- 
ered. Obviously, if the whole object is missing, it can not be 
recovered. 



4.3 Sparse representation of weak lensing mass 
maps 

Representing the image to be inpainted in an appropriate 
sparsifying dictionary is the main issue. The better the 
dictionary, the better the inpainting quality is. We are 
interested in a large and overcomplete dictionary that can 
be also built by the union of several sub-dictionaries, each 
of which must be particularly suitable for describing a 
certain feature of a structured signal. For computational 
cost considerations, we consider only (sub-) dictionaries 
associated to fast operators. Finding a sparse representation 
for weak lensing analysis is challenging. We want to describe 
well all the features contained in the data. The weak lensing 
signal is composed of clumpy structures such as clusters and 
filamentary structures (see Fig. [T]). The weak lensing mass 
maps thus exhibit both isotropic and anisotropic features. 
The basis that best represent isotropic objects are not the 
same as those that best represent anisotropic ones. We have 
consequently investigated a number of sub-dictionaries and 
various combinations of sub-dictionaries : 
- isotropic sub-dictionaries such as the "a trous" wavelet 
representation 
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Figure 6. Mean power spectrum error as a function of the max- 
imum number of iteration used by our inpainting method with 
the CFHTLS mask. 



Figure 5. Non-hnear approximation error I2 as a function of the 
percentage of coefficients used for the reconstruction, obtained 
with (i) the "a trous" Wavelet Transform (black), (ii) the local 
Discrete Cosine Transform (with a blocksize of 256 pixels) (red), 
(iii) the bi-orthognal Wavelet Transform (blue) and (iv) the "a 
trous" Wavelet Transform + the local DCT (blocksize = 256 pix- 
els) (green). The better representation for weak lensing data is 
obtained with a local DCT. 



- slightly anisotropic sub-dictionaries such as the bi- 
orthognal wavelet representation 

- highly anisotropic sub-dictionaries such as the curvelet 
representation 

- texture sub- dictionaries such as the Discrete Cosine 
Transform 

A simple way to test the sparsity of a representation 
consist of estimating the non-linear approximation error I2 
from complete data. It means, the error obtained by keeping 
only the N largest coefficients in the inverse reconstruction. 
Fig. [5] shows the reconstruction error I2 as a function of N. 

We expected wavelets to be the best dictionary because 
they are good to represent isotropic structures like clusters 
but surprisingly the better representation for weak lensing 
data is obtained with the DCT. More sophisticated repre- 
sentations recover clusters well, but neglect the weak lensing 
texture. Even combinations of DCT with others dictionar- 
ies (isotropic or not) is less competitive. We therefore chose 
DCT in the rest of our analysis. 



1. Set the maximum number of iterations Imax, the solution 

= 0, the residual to Y, the maximum threshold Xmax = 
max(|Q; = ^^Y\), the minimum threshold Xmin = 0- 

2. Set n to 0, An = Xmax- Iterate: 

3. U = X'^ + MR'' and R'' = {Y - X^). 

4. Forward transform of U: a = ^^U. 

5. Determination of the threshold level 

Xn = Fi^Tl, Xmax-, Xmin) ■ 

6. Hard-threshold the coefficient a using Xn : a = Sx^{a}. 

7. Reconstruct U from thresholded a and = <^q;. 

8. n = n and if n < Imax, return to Step 3. 



is the DCT operator. The way the threshold is de- 
creased at each step is important. It is a trade-off between 
the speed of the algorithm and its quality. The function F 
fixes the decreasing of the threshold. A linear decrease cor- 



responds to F(n, Xmax, Xmin) = \7nax Im,a^-1 

In practice, we use a faster decreasing law defined 

by. F {tI ^ Xmax 1 ^min) — ^min H~ (^^max '^m2n)(l- 

erf (2. Sn /Imax))- 

Constraints can also be added to the solution. For in- 
stance, we can, at each iteration, enforce the variance of the 
solution X^^^ to be equal inside and outside the masked 
region. We found that it improves the solution. 

The only parameter is the number of iterations Imax- 
In order to see the impact of this parameter, we made the 
following experiment : we estimated the mean power spec- 
trum error < Ep^{Imax) > for different values of Imax, see 
Fig.[6l The mean power spectrum error is defined as follows: 



4.4 Algorithm 

Our final dictionary $ being chosen, we want to minimize 
the number of non-zero coefficients (see eq. [9|). We use an 
i terati ve algorithm for image inpainting as in lElad et al.l 
that we describe below. This algorithm needs as in- 
puts the incomplete image Y and the binary mask M. The 
algorithm that we implement is: 



, (11) 



where k'T stands for the m*^ inpainted map with 
Imax iterations, Nq is the number of bins in the power spec- 
trum and Nm is the number of maps over which we estimate 
the mean power spectrum. 

It is clear that the error on the power spectrum de- 
creases and reaches a plateau for Imax > 100. We thus set 
the number of iterations to 100. 
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4.5 Handling noise 



If the data contains noise, it is straightforward to take it into 
account. Indeed, thresholding techniques are very robust to 
the noise since they are even used to remove it. In the pre- 
ceding algorithm, a filtered inpainted image can directly be 
obtained by setting the final threshold Xmin to ra instead 
of 0, where a is the noise standard deviation and r is a con- 
stant generally choosen between 3 and 5. However, even if 
the data is noisy, we may not want to perform denoising 
because it could introduce a bias in a further analysis such 
as power spectrum analysis. In this case, the final thresh- 
old Xmin should be kept at and the inpainting will try 
to reproduce also the noise texture (as if it were a real sig- 
nal). Obviously, it won't recover the "true" noise that will 
be observed if there was no missing data, but the statistical 
properties should be similar to that in the rest of the image. 
As we will see in the following, our experiments confirm this 
assertion. 



4.6 Experiments 

We present here several experiments to show how we have 
lowered the impact of the mask by applying the MCA- 
inpainting algorithm described above. The MCA-inpainting 
algorithm was applied with 100 iterations and using the 
DOT representation over the set of 100 incomplete mass 
maps, as described previously. The case of noisy mass maps 
has not been considered in the following experiments, it will 
be studied in ^5.41 



Incomplete mass map interpolation by MCA-inpainting 
algorithm 

The first experiment was conducted on two different simu- 
lated weak lensing mass maps masked by two typical mask 
patterns (see Fig. [T]). The upper panels show the simulated 
mass maps (see ^2.3p , the middle panels show the simulated 
mass map masked by the CFHTLS mask (on the left) and 
the Subaru mask (on the right) (see ^2.4|) . The results of 
the MCA-inpainting method using the DCT decomposition 
is shown in the lower panels allowing a first visual estimation 
of the quality of the proposed algorithm. We note that the 
gaps are no longer distinguishable by eye in the inpainted 
map. 



Probability Density Function comparison 

This second experiment was conducted on the weak lensing 
mass maps represented Fig|7] on the left. For these maps, 
we have computed the Probability Density Function (PDF) 
in order to compare their statistical distributions (see Fig. 
[8|). The PDF of the complete mass map is plotted as a solid 
black line, the PDF of the incomplete mass map as a solid 
blue line and the PDF of the inpainted mass map as a solid 
red line. The blue vertical line corresponds to the pixels 
masked out in incomplete mass maps. The strength of the 
PDF is that it provides a visual estimation of some statistics 
like the mean, the standard deviation, the skewness and the 
kurtosis. Thus, it provides a way to directly quantify the 
quality of the inpainting method. A visual comparison shows 



Complete mass map 
Inomplete mass map 
Inpainted mass map 




Figure 8. Probability Density Function estimated from the 3 
maps on the left of the Fig|7l : PDF of the complete simulated 
mass map (in black), PDF of the incomplete mass map (in blue) 
and the PDF of the inpainted mass map (in red). 



a striking similarity between the inpainted distribution and 
the original one. 



Power spectrum estimation 

This experiment was conducted over 100 simulated incom- 
plete mass maps. For these maps, for both the complete 
maps (i.e. no gaps) and the inpainted masked maps, we have 
computed the mean power spectrum: 



< P. >-- 



(12) 



where m is the number of simulations, the empirical stan- 
dard deviation also called (the square root of) sample vari- 



^^(P.^-<P.>)^. 



and the relative power spectrum error Ep : 



I\rr 



Pk"^ 



Pkr 



(13) 



(14) 



This experiment was done for the two kinds of mask 
(CFHTLS and Subaru). 

Fig.[9]shows the results for the CFHTLS mask. The left 
panel shows four curves: the two upper curves (almost su- 
perimposed) correspond to the mean power spectrum com- 
puted from i) the complete simulated weak lensing mass 
maps (black - continuous line) and ii) the inpainted masked 
maps (red - dashed line) . The two lower curves are the empir- 
ical standard deviation for the complete maps (black - con- 
tinuous line) and the inpainted masked maps (red - dashed 
hne). Fig. [9] right shows the relative power spectrum error, 
i.e. the normalized difference between the two upper curves 
of the left pannel. Fig. [TOl shows the same plots for the Sub- 
aru mask. 

We can see that the maximum discrepancy is obtained 
for / > 10000 with Subaru mask where the relative power 
spectrum error is about 3%. 
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Figure 7. Upper panels, simulated weak lensing mass map, middle panels, simulated mass map with the mask pattern of CFHTLS data 
on Dl field (left) and with the mask pattern of Subaru data in the same field (right), lower panels, inpainted mass map. The region 
shown is 1° x 1°. 



Computing time: 2PCF versus the inpainting-power 
spectrum method 

The two point-correlation function estimator is based on the 
notion of pair counting (see §[3T]). It is not biased by missing 



data, but its computational time is long. On our simulated 
field that covers a region of 1.975° x 1.975°, 8 hours are 
needed to process the two point correlation function in all 
the field with bins having the pixel size on a 2.5 GHz proces- 
sor PC-linux using C++. The proposed algorithm aims at 
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Figure 9. Power spectrum recovery from convergence maps for CFHTLS mask: left, the two upper curves (almost superposed) correspond 
to the mean power spectrum computed from i) the complete simulated weak lensing mass maps (black - continuous line) and ii) the 
inpainted masked maps (red - dashed line), and the two lower curves are the empirical standard deviation for the complete maps (black 
- continuous line) and the inpainted masked maps (red - dashed line). Right, relative power spectrum error, i.e. the normalized difference 
between two upper curves of the left pannel. 




Figure 10. Power spectrum recovery from convergence maps for the Subaru mask. 



lowering the impact of masked stars in the field while keeping 
fast calculation. The time to compute the power spectrum 
including the MCA-inpainting method in the same field still 
using the same 2.5 GHz processor PC-linux and C++ lan- 
guage is only 4 minutes. This is 120 times faster than the 
two-point correlation function. It only requires O(A^logA^) 
operations compared to the two-point correlation estimation 
that requires 0{N'^) operations. 



5 RECONSTRUCTION OF WEAK LENSING 
MASS MAPS FROM INCOMPLETE SHEAR 
MAPS 

In previous sections, we have investigated the impact of 
masking out the convergence field which is a good first ap- 
proximation. However, in real data, galaxy images are used 
to measure the shear field. This shear field can then be con- 
verted into a convergence field (see eg. 1151 and eg. I17|) . The 
mask for saturated stars is therefore applied to the shear 
field (i.e. the initial data), and we want to reconstruct the 
dark matter mass map k from the incomplete shear field 7. 

5.1 The inverse problem 

In weak lensing surveys, the shear ji{0) with i = 1,2 is 
derived from the shapes of galaxies at positions in the 



image. The shear field ji{0) can be written in terms of the 
lensing potential ip{0) as (see eg. iBartelmann and Schneider! 
(2001)): 

71 = ^ (^1 - ^2) ip 

72 = did2ip, (15) 

where the partial derivatives di are with respect to Oi. 

The projected mass distribution is given by the effective 
convergence k that integrates the weak lensing effect along 
the path taken by the light. This effective convergence can 
be written using the Born approximation of small scattering 
as (see eg. Bartelmann and Schneider (2001)): 

- 3Hin^ r Mw')Mw - w') s{Mw')d,w') ^ , 

Ke[ty) = — —5 — / . . aw , 

Jo fk[w) a{w') 

(16) 

where fk{w) is the angular diameter distance to the co- 
moving radius w, Hq is the Hubble constant, Qrn is the den- 
sity of matter, c is the speed of light and a the expansion 
scale parameter, 6 is the Dirac distribution. 

The projected mass distribution k,{0) can also be ex- 
pressed in terms of the lensing potential as : 

K = i (a? + di) i>. (17) 

The weak lensing mass inversion problem consists of re- 
constructing the projected (normalized) mass distribution 
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hi{0) from the incomplete measured shear field ji{0) by in- 
verting equations [15] and [171 This is an ill posed problem 
that need to be regularized. 



5.2 The sparse solution 

By taking the Fourier transform of equations [15] and [TT] we 

have 



(18) 



where the hat symbol denotes Fourier transforms. We define 
= -\- kl and 



A(k) 
A(k) 



fcf -fcj 

2fclfc2 
fc2 ' 



(19) 

when 



with Pi(/ci,/c2) = when ki — k^^ and P2{ki,k2) 
/ci = or /c2 = 0. 

We can easily derive an estimation of the mass map by 
inverse filtering, the least-squares estimator k of the conver- 
gence K IS e.g. (|Starck et al.ll2006l ): 



= Pi * 71 + P2 * 72- 



(20) 



We have ji = Pi^K, where * denotes convolution. When 
the data are not complete, we have: 



7, = M(P. */€), i = l,2. 



(21) 



To treat masks applied to shear field, the dictionary $ 
is unchanged because the DCT remains the best represen- 
tation for the data, but we now want to minimize: 

min ||$^A>:||o subject to ^ || ji-M(Pi^K) a. (22) 

i 

Thus, similarly to eq. [9] we can obtain the mass map k 
from shear maps ji using the following iterative algorithm. 



5.3 Algorithm 



1. Set the maximum number of iterations Imax, the solution 
= 0, the residual = Pi * 7^^^ + P2 * 7^^* see eq.[20]and 

^obs ^ (^o6s^^obs)_ ^he maximum threshold Xmax = max(|a = 
(f)^Y\), the minimum threshold Xmin = 0- 

2. Set n to 0, An = Xmax- Iterate: 

3. t/ = + MP^(7°^^) and 

P^(7^^^) = Pi * (jf^ - Pi * At^) + P2 * (jf^ - P2 * i^"") 

4. Forward transform of U: a = ^^U. 

5. Determination of the threshold level 

An = P(n,A 

max 5 Xmin ) ■ 

6. Hard-threshold the coefficient a using An : o; = Sx^{a}. 

7. Reconstruct U from the thresholded a and = ^a. 

8. n = n + 1 and if n < Imax, return to Step 3. 



5.4 Experiments 

One of the central goals of weak lensing analysis is the mea- 
surement of cosmological parameters. To constrain the large- 
scale structures, the power spectrum P^^ of the convergence 
K and thus its two-point correlation function C/^^, contains 
all the information about the primordial fluctuations. To 
characterize the non-gaussianity at small scales due to the 
growth of structures, higher-order statistics like three-point 
statistics have to be used. To lower the impact of the mask, 
we have applied the previous algorithm with 100 iterations 
and using always one single representation, the DCT. Then 
we have conducted several experiments to estimate the qual- 
ity of the method. 



Dark matter mass map reconstruction from incomplete 
shear maps 

As in the previous section, we have applied the MCA- 
inversion method on two different simulated weak lensing 
mass maps whose shear field have been masked by the two 
typical mask patterns. Fig. [11] top panels, shows the simu- 
lated mass map masked by the CFHTLS mask (on the left) 
and by the Subaru mask (on the right). The results apply- 
ing the MCA-inpainting method described above is shown 
in the bottom panels. The gaps are again undistinguishable 
by eye in both cases. 

Power spectrum estimation of the convergence k, from 
incomplete shear maps 

As in the previous section, this experiment was done over 
100 incomplete shear maps for the two kinds of mask 
(CFHTLS and Subaru). 

Fig. [12] shows the results for the CFHTLS mask. The 
left panel shows four curves: the two upper curves (almost 
superposed) correspond to the mean power spectrum com- 
puted from i) the complete simulated weak lensing mass 
maps (black - continuous line) and ii) the inpainted recon- 
structed maps from masked shear maps (red - dashed line). 
The two lower curves are the empirical standard deviation 
estimated from the complete maps (black - continuous line) 
and the inpainted reconstructed maps from masked shear 
maps (red - dashed line). Fig. [12] right shows the relative 
power spectrum error, i.e. the normalized difference between 
the two upper curves of the left pannel. And the blue - 
dashed line is the root mean square of the sample variance 
that comes from the finite size of the field. Fig. ll3lshows the 
same plots for the Subaru mask. 

We can see that the maximum discrepancy is obtained 
in the /-range of [2000, 7000] with CFHTLS mask where 
the relative power spectrum error is about 1% while is 
about 0.3% with Subaru mask. The error introduced by our 
method is consistent with the sample variance. 



is the DCT operator. The residual Rn is estimated 
from shear maps 7^^* and 72^*. Consequently, we need to 
use two FFTs at each iteration n to compute the mass map 
K from the shear fields (eq. I20|) and the shear fields from 
the mass map k fea. I18|) . F follows the same decreasing law 
described in ^4.41 



Case of noisy incomplete shear maps 

As discussed in section [43] if the data contains some noise, 
it is straightforward to take it into account using a simi- 
lar processing. The weak lensing data are noisy because the 
observed shear ji is obtained by averaging over a finite num- 
ber of galaxies. Another experiment has been conducted still 
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Figure 11. Upper panels, simulated weak lensing mass maps with the mask pattern of CFHTLS data on Dl field (left) and with the 
mask pattern of Subaru data in the same field (right) whose masks for saturated stars are applied to the shear field. Lower panels, 
inpainted mass maps. The region shown is 1° x 1°. 



using 100 simulated weak lensing mass maps with noise to 
simulate space observations (100 galaxies/arcmin^ and with 
the shear error per galaxy given by = 0.3). The two kinds 
of mask (CFHTLS and Subaru) have been also tested. 

Fig. [H shows the results for the CFHTLS mask. Left 
panel, the two upper curves (almost superposed) correspond 
to the mean power spectrum computed from i) the complete 
simulated noisy mass maps (black - continuous line) and ii) 
the inpainted maps from masked noisy shear maps (red - 
dashed line) . The noise power spectrum represented in green 
dashed line modify the shape of the noiseless weak lensing 
power spectrum plotted in blue dashed line. The two lower 
curves are the empirical standard deviation for the complete 
noisy mass maps (black - continuous line) and the inpainted 
maps from masked noisy shear maps (red - dashed line). 
Fig. [TJ] right shows the relative power spectrum error, i.e. 
the normalized difference between two upper curves of the 
left pannel. Fig. [15] shows the same plots for the Subaru 
mask. 

We can see that the maximum discrepancy is obtained 
with CFHTLS mask in the /-range of [25000, 40000] where 
the relative power spectrum error is about 1% and is only 



0.5% with Subaru mask. The error is not amplified by the 
noise. 

Equilateral hispectrum estimation of the convergence k 
from incomplete shear maps 

We consider now the equilateral bispectrum that we have 
introduced ^3.21 We defined the mean equilateral bispectrum 
as follows : 

<B:'>=-)-^B:L. (23) 

m 

And the relative equilateral bispectrum error E^eq is given 
by: 

The experiment was also done for the two kinds of mask 
(CFHTLS and Subaru). Fig. 1161 shows the results for the 
CFHTLS mask. The left panel shows three curves (whose 
two almost superposed) that correspond to the mean equi- 
lateral bispectrum computed from i) the complete simulated 
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Figure 12. Power spectrum recovery from shear maps for CFHTLS mask: left, the two upper curves (almost superposed) correspond 
to the mean power spectrum computed from i) the complete simulated weak lensing mass maps (black - continuous line) and ii) the 
inpainted reconstructed maps from masked shear maps (red - dashed line), and the two lower curves are the empirical standard deviation 
for the complete maps (black - continuous line) and the inpainted reconstructed maps from masked shear maps (red - dashed line). 
Right, relative power spectrum error, i.e. the normalized difference between two upper curves of the left pannel. The blue - dashed line 
represents the empirical standard deviation (cosmic variance) estimated from the complete mass maps. 




Figure 13. Power spectrum recovery from shear maps for the Subaru mask. 



mass maps (black - continuous line) and ii) the inpainted 
maps from masked shear maps (red - dashed line) and iii) 
the incomplete simulated shear maps (green - continuous 
line). Fig. [16] right shows the relative equilateral bispectrum 
error in logarithmic bins, i.e. the normalized difference be- 
tween the curves of the left pannel. Fig. 1171 shows the same 
plots for the Subaru mask. 

We defined the mean equilateral bispectrum as follows 

<B7>=^Y.B7^- (25) 



And the relative equilateral bispectrum error E^eq is given 
by: 




The maximum discrepancy is obtained with the 
CFHTLS mask where the relative bispectrum error is about 
3% while is about 1% with Subaru mask, which remains con- 
sistent with the sample variance in blue - dashed line on the 
right. This result is satisfactory because no constraint on the 
solution has been used to improve the estimation quality of 
the bispectrum. 



Computing time: two-point correlation versus 
power spectrum using the iterative reconstruction 

As in the previous section, we have compared the comput- 
ing time of the two-point correlation function applied on the 
incomplete shear maps with the power spectrum applied on 
inpainted mass map obtained from incomplete shear maps. 
The time to process the MCA-inpainting starting from shear 
maps followed by the power spectrum is 6 minutes still us- 
ing a 2.5 GHz PC-linux processor. It is a bit longer than 
starting from convergence maps because it requires a con- 
version from shear field to convergence field and vice versa 
at each iteration in the MCA-inpainting and also because 
the algorithm is this time written in IDL language. It still 
remains 80 times faster than the two-point correlation func- 
tion. Furthermore, the reconstructed map can also be used 
to compute higher-order statistics. 

Computing time : three-point correlation versus 
inpainting reconstruction followed by a bispectrum 

Here, we compare the time needed to compute on the one 
hand the three-point correlation function applied on the in- 
complete mass maps and on the other hand, the bispectrum 
applied on inpainted masked shear maps. Our three point- 
correlation function estimator is based on the averaging of 
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Figure 14. Power spectrum recovery from noisy shear maps for CFHTLS mask: left, the two upper curves (almost superimposed) 
correspond to the mean power spectrum computed from i) the complete simulated noisy mass maps (black - continuous line) and ii) 
the inpainted reconstructed maps from masked noisy shear maps (red - dashed line). The mean power spectrum of the noise is plotted 
as a dashed green line and the mean power spectrum of the noiseless mass maps is plotted as a dashed blue line. The pink dashed 
line represents the simulation shot noise. The two lower curves are the empirical standard deviation for the complete noisy mass maps 
(black - continuous line) and the inpainted reconstructed maps from masked noisy shear maps (red - dashed line). Right, relative power 
spectrum error, i.e. the normalized difference between two upper curves of the left pannel. 
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Figure 15. Power spectrum recovery from noisy shear maps for the Subaru mask. 



the power of the equilateral triangles on the direct space and 
it is not biased by missing data. But its computational time 
is very long. Even using the equilateral information, more 
than 8 hours are needed to process the three-point corre- 
lation function in a field of 4 square degree on a 2.5 GHz 
processor PC-linux. 

The time to compute the equilateral bispectrum includ- 
ing the MCA-inpainting method all written in IDL, in the 
same field still using the same 2.5 GHz processor PC-linux is 
only 6 minutes. This is 80 times faster than the three-point 
correlation function. And it only requires O(A^logA^) op- 
erations compared to the three-point correlation estimation 
that requires 0{N'^) operations. 



6 CONCLUSION 

This paper addresses the problem of statistical analysis of 
Weak Lensing surveys in the case of incomplete data. 

We have presented a method for image interpolation 
across masked regions and its applications to weak lens- 
ing data analysis. The proposed inpainting approa ch relies 
strongly on the ideas of MCA (jStarck et al.ll2004l l and on 
the sparsity of the weak lensing signal in a given dictionary. 
The proposed reconstruction algorithm is based on a decom- 



position in a basis of cosines that turned out to be the best 
representation for weak lensing data. This recent tool can 
be applied to many other interesting applications and it is 
alrea dy used in CMB da ta analysis to fill in the galactic re- 
gion ([Abrial et al.|[2007l ). The algorithm has been extended 
to the weak lensing inversion problem with realistic masks. 
With this extension, the inpainting method provides a solu- 
tion for the problem of estimation of the convergence mass 
map from incomplete shear maps. 

We have proposed to use our fast 0{N log N) inpainting 
algorithm to lower the impact of missing data on statistics 
estimations. We have shown that our inpainting method en- 
ables us to use the power spectrum in future large weak 
lensing surveys by filling in the gaps in weak lensing mass 
maps in the presence of noise. We have shown that our in- 
painting method enables to reach an accuracy of about 1% 
with the CFHTLS mask and about 0.3% with Subaru mask 
for the power spectrum 

We have shown that our inpainting technique can also 
be applied to higher-order statistics. In particular, we have 
presented a fast and accurate method to calculate the equi- 
lateral bispectrum using a polar FFT and we have shown 
that our inpainting method enables to reach an accuracy 
of about 3% with the CFHTLS mask and about 1% with 
Subaru mask for the bispectrum. 
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Figure 16. Bispectrum recovery for CFHTLS mask: left, the two upper curves (almost superposed) correspond to the mean equilateral 
bispectrum computed from i) the complete simulated weak lensing mass maps (black - continuous line) and ii) the inpainted maps from 
masked shear maps (red - dashed line), and the lower curve correspond to the mean equilateral bispectrum computed from the incomplete 
shear maps (green - continuous line). Right, relative power spectrum error, i.e. the normalized difference between the curves of the left 
panel in logarithmic bins. The blue - dashed line represents the empirical standard deviation estimated from the complete mass maps, 
previously on black - continuous line. 
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Figure 17. Bispectrum recovery for the Subaru mask. 



It w ould be interesting to extend the MRLENS filtering 
package (|Starck et al.ll2006 l in order to build a filtered mass 
map from incomplete shear maps by applying the inpainting 
technique. 

In future work, this inpainting technique can be ex- 
tended to compute other non-gaussian statistics : higher- 
order statistics, peak finding, etc... taking advantage of the 
map recovery. 



Software 

The software related to this paper, FASTLens, and its full 
documentation is available from the following link: 

\protect\vrule widthOpt\protect\href {http : // j 
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